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We study the non equilibrium dynamics in the fermionic Hubbard model after a sudden change 
of the interaction strength. To this scope, we introduce a time dependent variational approach in 
the spirit of the Gutzwiller ansatz. At the saddle-point approximation, we find at half filling a 
sharp transition between two different regimes of small and large coherent oscillations, separated 
by a critical line of quenches where the system is found to relax. Any finite doping washes out 
the transition, leaving aside just a sharp crossover. In order to investigate the role of quantum 
fluctuations, we map the model onto an auxiliary Quantum Ising Model in a transverse field coupled 
to free fermionic quasiparticles. Remarkably, the Gutzwiller approximation turns out to correspond 
to the mean field decoupling of this model in the limit of infinite coordination lattices. The advantage 
is that we can go beyond mean field and include gaussian fluctuations around the non equilibrium 
mean field dynamics. Unlike at equilibrium, we find that quantum fluctuations become massless 
and eventually unstable before the mean field dynamical critical line, which suggests they could 
even alter qualitatively the mean field scenario. 

PACS numbers: 71.10.Fd, 05.30.Fk, 05.70.Ln 



I. INTRODUCTION 

Recent years have seen an enormous progress in 
preparing, controlling and probing ultra cold atomic 
gases loaded in optical lattices^. Their high degree of 
tunability allows to change in time the microscopic pa- 
rameters controlling interactions among atoms and to 
measure the resulting quantum evolution. At the same 
time their excellent isolation from the environment makes 
those systems particularly well suited to address ques- 
tions related to non equilibrium phenomena in isolated 
many body quantum systems. These major achievements 
triggered a huge interest on time dependent phenomena 
in condensed matter systems. In this respect, the re- 
cent experimental realization of a fermionic Mott insu- 
lator^i 3 - opened the way to investigate out-of-equilibrium 
phenomena in strongly correlated fermionic systems^. 

From a more theoretical perspective these experiments 
offer the chance to probe strongly correlated systems 
in a completely novel regime. Indeed, when driven out 
of equilibrium, interacting quantum systems can display 
peculiar dynamical behaviors or even be trapped into 
metastable configurations that differ completely from 
their equilibrium counterpart 5,6 . Although actual exper- 
iments are always performed by tuning parameters at a 
finite rate, an useful idealization consists in a so called 
quantum quenchZ. Here the system is firstly prepared in 
the many-body ground state \^i) of some initial Hamil- 
tonian Hi which is then suddenly changed to Hf ^ Hi, 
for example by globally switching on or off some coupling 
constants. As a consequence of this instantaneous change 
the initial state \^i) turns to be an highly excited state 



of the final Hamiltonian. Naturally, many non trivial 
questions arise concerning the real-time evolution after 
the quantum quench. The interest on these classes of 
non equilibrium problems relies both on the dynamics 
itself*^ as well as on the long-time properties where the 
question of thermalization or its lack of is still highly 
debated^Sr— This issue is not only of fundamental the- 
oretical interest but also of great practical relevance for 
establishing whether and to what extent experiments on 
cold atoms could reproduce equilibrium phase diagrams 
of model hamiltonians. 

The literature on quantum quenches in interacting 
bosonic and fermionic systems is by now very broad, see 
for example the recent topical reviews^ 3 - - — For what con- 
cerns strongly correlated electrons in more than one di- 
mension, the subject is still largely unexplored and pro- 
gresses have been done only very recently. The single 
band Hubbard modeli^r— represent one of the simplest 
yet non trivial models encoding the physics of strong 
correlations, namely the competition between electronic 
wave function derealization due to hopping t and charge 
localization due to large Coulomb repulsion U. Its Hamil- 
tonian reads 

W (*) = - 53 53 tRR ' C R<7 C R'a + U{t)Y, n ^t ™IU • 
<T (R,R') R 

Despite the everlasting interest on its groundstate prop- 
erties, theoretical investigations on the non equilibrium 
dynamics of this paradigmatic strongly correlated model 
have been started only very recently. The dynamics of 
Fermi system after a sudden switch-on of the Hubbard in- 
teraction has been studied firstly in Refs. I20II21I using the 
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flow-equation approach and then in Refs . [22I I 23I using Non 
equilibrium Dynamical Mean Field Theory (DMFT) . Re- 
sults suggest the existence of two different regimes in the 
real-time dynamics, depending on the final interaction 
strength Uf. At weak coupling^ the systems is trapped 
at long-times into a quasi-stationary regime which looks 
as a zero temperature Fermi Liquid from the energetic 
point of view but features a non thermal distribution 
function in which correlations are more effective than in 
equilibrium. This pre-thermalization phenomenon has 
been confirmed by DMFT results* 2 ^ which further in- 
dicate a true dynamical transition above a critical U f c 
towards another regime with pronounced oscillations in 
the dynamics of physical quantities. This picture has 
been recently confirmed by means of a simple and flex- 
ible approximation scheme based on a proper extension 
of the Gutzwiller variational method 24 . Results for the 
time depedent mean field theory show at half filling, a 
sharp transition between two different regimes of small 
and large coherent oscillations, separated by a critical 
line of quenches where the system finds a fast way to re- 
lax. Away from particle hole symmetry the transition is 
washed out, leaving a sharp crossover visible in the dy- 
namics and in the long-time averages of physical quanti- 
ties. 

The aim of the present work is twofold. From one 
side, we present details on the time dependent Gutzwiller 
method for fermions and discuss its application to the 
problem of an interaction quench in the single band Hub- 
bard model. Secondly, we discuss the role of quantum 
fluctuations on top of the Gutzwiller dynamics. In order 
to do that we formulate the original Hubbard model in 
terms of an auxilary Quantum Ising Model in a trans- 
verse field coupled to free fermionic quasiparticle. Such 
a Z2 slave spin theory, introduced in Refs. [25ll26l for the 
equilibrium problem, allows us to study the effect of small 
quantum fluctuations, both in equilibrium as well as for 
the non equilibrium dynamics. We notice that the role of 
quantum fluctuations on this mean field dynamical tran- 
sition is of broader theoretical interest, as recent investi- 
gations have shown the very same phenomenon occurs in 
other models of interacting quantum field theorie a 27 ' 28 . 

The paper is organized as follows. In the first part 
we introduce the time dependent variational method we 
have devised to describe non equilibrium dynamics in 
correlated electrons systems. Section |TT] is devoted to 
a general formulation while section Mil to the study of 
quantum quenches in the single band fermionic Hubbard 
Model. In the second part of the paper we broaden the 
perspective and formulate the Hubbard model in terms of 
auxiliary Quantum Ising Model coupled to free fermionic 
quasiparticles. In section ITVl we show how the mapping 
works and how to recover the Gutzwiller results. Sec- 
tion [IVD] is devoted to the role of quantum fluctuations. 
Finally section [V] is for conclusions. 



II. A GENERAL FORMULATION 

We assume a system of interacting electrons that is 
initially in a state with many-body wavefunction \^o). 
For times t > 0, |\&o) is let evolve with the Hamiltonian 
H, which could even be explicitly time-dependent. We 
shall assume that short range correlations are strong ei- 
ther in the initial wavefunction, or in H, or in both. The 
goal is calculating average values of operators during the 
time evolution. Because of interaction, a rigorous calcu- 
lation is unfeasible, so that an approximation scheme is 
practically mandatory. Our choice will be to use a proper 
extension of the Gutzwiller wavefunction and approxima- 
tion, which is known to be quite effective at equilibrium 
when strong short-range correlations are involved. 

We start by defining a class of many-body wavefunc- 
tions of the form 

!*(*)> = l[e- iS ^V R (t)m)) 

R 

= P(i)|$(t)>, (2) 

where |$(f)) are time-dependent variational wavefunc- 
tions for which Wick's theorem holds, hence Slater de- 
terminants or BCS wavefunctions, while Vji(t) and 6>rq, 
are hermitian operators that act on the Hilbert space at 
site i and depend on the variables Ar,q,(£) and </>R a (i): 

Pn(t) = W*)0Ra> (3) 

Ra 

-^-e- i5 * = - 8 Ra e- ,& , (4) 

0(f>R a 

where Or q can be any local hermitian operator. It fol- 
lows that the average value of Or q 

Or q = (*(t)|0 RQ |*(t)), (5) 

is a functional of all the variational parameters. We shall 
assume that it is possible to invert © and express the 
parameters Ar q as functionals of all the Or'^, 0r'/j as 
well as of the parameters that define 

Since ^(t)) spans a sub-class of all possible many- 
body wavefunctions, in general it does not solve the 
Schrcedinger equation but can be chosen to be as close 
as possible to a true solution. This amounts to search for 
the saddle point of the functional 

S[*t,tf] = J dt(*{t)\id t -H\*(t)), (6) 

with of the form as in Eq. fl2J. The Gutzwiller ap- 

proximation gives a prescription for calculating S, which 
is exact in infinite coordination lattices ; 29 ' 30 although it 
is believed to provide reasonable results also when the 
coordination is finite. We impose that 

m)\-pi(t)mt)) = 1, (7) 

p R (t)c Ra |<i>(*)) = <$(t)|c RQ !<&(*)>, (8) 
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where Cr q is any bilinear form of the single-fermion op- 
erators at site R, c Ra and c Ra with a the spin/orbital 
index. 

Within the Gutzwiller approximation and provided 
Eqs. and (jHJ) hold, the average value of any local op- 
erator Orq, is assumed to be^ 

Ra = (*(t)\0 Ra \V(t)) = (9) 



- (*(t)| W) e lSR(t) Ra e" iS *« V R (t) |*(t)), 



which can be easily computed by the Wick's theorem. 
Seemingly, given two local operators, Orq, and Or'/j at 
different sites R ^ R', the following expression is as- 
sumed 



(*(t)| Or q O wp |*(t)) - ($(*)| Pn(t) e l5R(t) Or q e- i5 ^( { ) Pn(t) pR,(t) e i5R ' W Or^ e" 



iS R ,(t) 



!*(*)>, (io) 



r 



which can be also readily evaluated. For consistency, one 
should keep only the leading terms in the limit of infinite 
coordination lattices^ For instance, if |$(t)) is a Slater 
determinant and Or„ = c Ra while 0R<f, = c R , b , then 

(* (*)l 4a °R> 6 !*(*)> = 

= EQR.«cQ R ', M ^(*)l4cCR' d l $ (*)): (11) 

where the matrix elements Qn^b are obtained by solving 
<*(t)| Pr(*) e lS - (t) 4 a e- <Sa(i) PrW c Rc |$(t)> 

= E 4* Crc !*(*)>• (i 2 ) 

b 

Within the Gutzwiller approximation one finds that 

^(*(^)|^ i *(^)) - J2 fa°> + *<*(*)! W*)>> (1 3 ) 



Rn 



so that 



5[* t ,*] = y <ft <j>RaO Ra -E[<f> 

Rq 



(14) 



£[fea,OR ai $] = (15) 

ft, = P\t)UP(t). (16) 



where 



The saddle point of 5 in Eq. (fT4| with respect to 0r o 
and Or q is readily obtained by imposing 



• _ dE 
■ dE 

° Ra — — ■ 

OPRa 



(17) 
(18) 



showing that these pairs of variables act like classical 
conjugate fields with Hamiltonian E. As far as 



is concerned, since it is either a Slater determinant or a 
BCS wavefunction, the variation with respect to it leads 
to similar equations as in the time-dependent Hartree- 
Fock approximation, 31 namely in general, non-linear sin- 
gle particle Schrcedinger equations. 

In conclusion, the variational principle applied to the 
Schrcedinger equation and combined with the Gutzwiller 
approximation amounts to solve a set of equations that 
is only slightly more complicated than the conventional 
time-dependent Hartree-Fock approximation, yet incom- 
parably simpler than solving the original Schrcedinger 
equation. We note that, in the above scheme, the 
Gutzwiller variational parameters Ar q in Eq. or 
better Or q in Eq. (jSJ), have their own dynamics be- 
cause of the presence of their conjugate fields </>r q . This 
marks the difference with the time-dependent variational 
scheme introduced by Seibold and Lorenzana 32 , where 
the time evolution of Ar q is only driven by the time evo- 
lution of the Slater determinant. We shall see that this 
difference may play an important role. 



III. QUANTUM QUENCHES IN THE 
HUBBARD MODEL 



We now turn to the problem of our interest and discuss 
the non equilibrium dynamics in the Hubbard model (JTJ 
using the time dependent variational scheme introduced 
above. This calculation allows to benchmark the method 
towards more reliable techniques, a compulsory step be- 
fore moving to more complicated situations where rigor- 
ous results are lacking. In particular we shall study the 
dynamics after a sudden change of the local interaction, 
starting from the zero-temperature variational ground 
state with U(t < Q) = f/j then quenching the interac- 
tion to U(t > 0) = Uf. Notice that since the initial state 
is described within the equilibrium Gutzwiller approx- 
imation, which provides a poor description of the Mott 
Insulator, we have to restrict our analysis to strongly cor- 
related yet metallic initial conditions, namely to Ui < U c 
where U c is the critical interaction strength for the Mott 
transition within the Gutzwiller approximation. More- 
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over, in what follows we shall completely disregard mag- 
netism, considering only paramagnetic and homogeneous 
wave functions. 



A. Time Dependent Gutzwiller Approximation 

We take T-L to be the single band Hubbard model ([1]) 
and assume a correlated time-dependent wave function 
of the form (|2|) with 



V R (t) = A R , n (t)P R , n , 

n=0 
2 

5 R (t) = fo,»(*)PR,n, 



(19) 



(20) 



n=0 



where Vn.n is the projector at site R onto configurations 
with n = 0, . . . , 2 electrons. Notice that equations (| 194 - 
|20| imply that 4>B,,n(t) plays the role of the conjugate 
variable of 



PR,n = (^(t)\V R:n \^(t)). 



(21) 



For non-magnetic wave functions, the renormalization pa- 
rameters in Eq. (|12j) do not depend on the spin index and 
read 



Qi = 



where 



is the average on-site occupancy. The two constraints 
Eqs. (0) and dHJ imply that the quantities Pr,„ in (|2"Tj) 
behave as genuine occupation probabilities with 

n 

If we set Pr :2 = Br then Pro = 1-?^r+-D r and Prs = 
tir — 2Z3r. We also assume that <^r,o = <^r,2 — '/'R while 
i^r^ = 0, so that the energy functional E becomes 

£[0r,Pr,$] = (V(t)\H\V(t)) =U f J2 D * + 

R 

+ E Q R QH/«*lR'(*) + ff-C#3) 
(RR'j 



where 



while Qr(<) defined in equation (|22j) reads 



n R - 2Z> 



r 



R 



n-R (1 - nR/2) 



(/Dr~ 



1 - «r e l 



(25) 



By the variational energy (|23|) we can readily obtain the 
equations of motion for the double occupancy Dr and 
its conjugate variable </)r using (|17I18I) . In addition, 
the dynamics of these variational parameters is further 
coupled to a time dependent Schroedinger equation for 
the Slater determinant. If this latter is initially homo- 
geneous, then translational symmetry is mantained dur- 
ing the time evolution, hence Qr(£) = Q(t) independent 
of R. Moreover, if the Slater determinant |$(i = 0)) 
is initially the Fermi sea, i.e. the lowest energy eigen- 
state of the hopping Haimiltonian, then its time evolu- 
tion caused by the time dependent hopping | Q{t)\ 2 £rr< 
becomes trivial 

\$(t))=exp(-iVe n ^dT\Q(T)\ 2 ^ |3(i)>, 

where e n is the average energy per site of the hopping 
Hamiltonian with electron density nona lattice with V 
sites. In other words, the matrix elements wrr< (t) in 
Eq. (j2~4")l are in this case time independent. 



B. Saddle-point equations 

In conclusion, within the Gutzwiller approximation 
and assuming a homogeneous and non-magnetic wave- 
function the classical Hamiltonian (l23l) for the single 
degree of freedom Dr = D and its conjugate variable 
</>r = (j) reads 



E[D,4>] = U f D(t)+e n Z(D, 



(26) 



where we remind that e n is the average hopping energy of 
a Fermi sea with density n = 1 — 6 while Z = \Q\ 2 is the 
effective quasiparticle weight, which reads from equation 



2 (n - 2D) 
n (2 - n) 



(27) 



VW+5 - Vd) + 4 cos 2 (bVD^D + S 



Z(D,, 



Notice that Z does not depends only from the double 
occupation D, as one would expect in equilibrium, but 
features a dependence from the phase <f> which is crucial 
in order to induce a non trivial dynamics. 

The classical equations of motion for this intcgrable 
system immediately follow from 
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S_ndZ_ 

2 dD 1 



(28) 
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FIG. 1: Sketch of the phase diagram in the Ui,uf plane for the 
quench dynamics of the single band Hubbard model within 
the Gutzwiller approximation at half-filling. Two different 
dynamical regimes corresponding to weak and strong coupling 
dynamics (A and B in the plot) are found depending whether 
the final interaction Uf lies above or below the critical quench 
line Uf c — 1+ 2 i ■ For quantum quenches along this line the 
dynamics features an exponential relaxation toward a steady 
state. 
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e_ndZ_ 
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(29) 



In the following we will use the MIT critical interaction at 
half-filling, U c = — 8e„ = i = — 8e, as the basic unit of en- 
ergy and define accordingly the dimensionless quantities 
Uf = Uf/U c and Ui = Ui/U c , as well a dimensionless time 
t = tU c . In addition we shall assume for simplicity a flat 
density of states so that e„ = n(2 — n)e — —n(2 — n)/8U c . 

The initial conditions for the classical dynamics 
(1291) read 



D{t = 0) = A , <f>(t = 0) = , 



(30) 



where Di is the equilibrium zero temperature double oc- 
cupancy for interaction m and doping 6 that can be eas- 
ily computed from an equilibrium Gutzwiller calculation, 
which is nothing but annihilating the right hand sides of 
Eqs. (|28|l and (|2"9"|) with interaction Ui instead of Uf. 

It is worth noticing that, apart from the trivial case in 
which Uf = Ui, the classical dynamics (| 2 8 1) - ([29"]) admits 
a non-trivial stationary solution D = and cos 2 6 = Uf, 
which is compatible with the initial conditions only at 
half- filling and Uf = Uf c = (1 + ttj)/2. It turns out 
that Uf c identifies a dynamical critical point that sepa- 
rates two different regimes similarly to a simple pendu- 
lum. When Uf < Uf c , 26{t) oscillates around the origin, 
while, for Uf > Uf c , it performs a cyclic motion around 
the whole circle. In order to characterize the different 
regimes, we focus on three physical quantities, the dou- 
ble occupancy D(t), the quasiparticle residue Z(t) and 
their period of oscillation, T. 



Before discussing in some detail the results of the clas- 
sical dynamics (I2"8l) - ([2U)) . it is useful to cast it into a closed 
first-order differential equation for one of the two conju- 
gate variables D or (p. Indeed the dynamics conserves 
the energy, namely 

n(2 — n) 

E(t)=u f D{t) ^ — -Z(t) = E , t>0 (31) 

8 

where Eg is the total energy soon after the quench, which 
reads 



E = u f Di 



n(2 - n) 



(32) 



with Zi the equilibrium zero temperature quasiparticle 
weight for interaction Ui and doping S. The simplest way 
to proceed is to eliminate 6 from Eq. (l3"TT) in favor of the 
double occupancy D{t). From Eq. (j27|) we obtain 



E - u f D + (n - 2D) ( y/17+6 - y/D) /4 



COS = — - 



(n - 2D) D (D + 5) 



(33) 

which can be inserted into (|29l) and leads, after some 
algebra, to the equation of motion 



D = ±^f{D). 



(34) 



Here T(D) can be thought as an effective potential con- 
trolling the dynamical behavior of D(t). We note that, 
since the problem is one dimensional, many properties of 
the solution (|34|) can be inferred directly from the knowl- 
edge of r(D), without explicitly solving the dynamics. In 
the next two sections we will discuss in detail the struc- 
ture of this solution, considering both the half filled and 
the doped case. 



C. Quench Dynamics at Half-Filling 

We start by considering half- filling, i.e. 6 = 0, and 
for simplicity we fix Uf > Ui, see figure [TJ As we al- 
ready anticipated, the dynamical behavior of the system 
changes drastically when the final value of the interaction 
Uf crosses the critical line Uf c = (1 + ui) /2. 

The existence of such a line of critical values clearly 
emerges from the structure of the effective potential 
r (D) and in particular from the behavior of its positive 
roots, which are the inversion points of the one dimen- 
sional motion (|3"4"|) . 

As one can see from figure [21 T(D) has three simple 
zeros, two of them being positive. It turns out that the 
equilibrium Gutzwiller solution Di is always one of the 
roots of the effective potential, for any Uf, see figure [2] 
(top panels). The remaining two, D±, depend strongly 
on Uf as we show in the bottom panel of figure [5] Since 
the one dimensional motion is constrained to the interval 
\D + , Di], where r(D) is positive, we expect to find peri- 
odic solution of the dynamics (|34[) . However, the proper- 
ties of this solution will largely depend on the behavior 
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FIG. 2: Top Panel: Effective potential V (D) for u f = 
0.2,0.3,0.4,0.5 (right) and u f = 0.6,0.8,0.9,1.0 (left). Bot- 
tom panel: Inversion points D+,D_ as a function of u/ at 
fixed Ui = for zero and finite doping. 



FIG. 3: Left Panel: Mean-field dynamics for quantum 
quenches to u/ below (top) and above (bottom) the critical 
line. Right Panel: Period of oscillation Tb as a function of w/ 
for m — 0.0, 0.25, 0.5, 0.75. Notice the log-singularity at u/ c . 



of D + as a function of Uf. As we see, D + first decreases 
linearly with Uf, vanishing at Uf c where it becomes de- 
generate with D_ (see figured]). Then for Uf > Uf c it 
starts increasing again, approaching Di in the infinite 
quench limit. It turns out that D + has a simple form, 
which reads 




■/2 



"f. 



(35) 



Two different dynamical behaviors are therefore ex- 
pected as a result of this peculiar dependence. In ad- 
dition, due to the degeneracy of simple roots occurring 
at Uf — Uf c we expect here a special trajectory, where 
relaxation to a steady state can exist. This qualitative 
picture is confirmed by the actual solution of the classi- 
cal dynamics (jMJ), whose results we are going to present, 
both for weak (uf < Uf c ) and strong (uf > Uf c ) quantum 
quenches. 



Weak Quenches: Uf < uf c 

For weak quantum quenches to Uf < Uf c , the dynamics 
of both double occupation D(t) and quasiparticle weight 
Z(t) shows coherent oscillations, see figure [3l which do 
not die out. The lack of relaxation toward a steady state 
is clearly an artifact of our semiclassical approach that 
does not account for quantum fluctuations. This is par- 
ticularly true for weak quenches starting from the gapless 
metallic phase, where fast damping of the oscillations is 
expected due to the available continuum of low-lying ex- 
citations. 



Although oversimplified, the dynamics in the weak 
quench limit contains some interesting features that are 
worth to discuss. In particular, we focus on the period 
T of the coherent oscillations as a function of the final 
interaction Uf. It is easy to see that T is given by 



T=2 



dD _4V2K(k) 



(36) 



where K (k) is the complete elliptic integral of the first 
kind with argument k 2 — Auf (uf — ui) jZi. As we show 
in the right panel of figure [3l upon increasing u / the 
period T grows eventually diverging logarithmically as 
the critical quench line Uf = Uf c is approached. This 
can be seen explicitly in Eq. (l36[) . Indeed, for Uf — ¥ Uf c , 
the argument of the complete elliptic integral approaches 
k = 1 



1 



{U fc -Uf)(l + 



2u fc D. 

Therefore, using the known asymptotic result K (k) 
log (4/Vl - k 2 ) we find 



(37) 



log 



1 



Uf c - Uf 



(38) 



Such a diverging time scale signals a sharp transition to 
a completely different dynamical regime for Uf > Uf c . 
Before moving to this strong coupling regime we briefly 
discuss the dynamics of the phase 4>{t) in the weak quench 
case, which can be easily obtained by eliminating the 
double occupation D(t) from the original system (|28H29I) . 
As shown in figure 21 in the present weak quench regime 
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FIG. 4: Dynamics of the phase in for quenches below and 
above the critical value uf c . Notice that for small quenches 
the phase oscillates around zero while for Uf > Uf c the dy- 
namics is no more bounded since the energy is sufficient to 
overcome the potential barrier. 



(v.f < Uf c ) the phase oscillates around the equilibrium 
fixed point = 0, with the same period T . As we are 
going to discuss in the next paragraph, it is just the phase 
which shows the most striking change in the dynamics as 
the critical line is crossed. 



Strong Quenches: it/ > u/ c 

As we anticipated, for quenches above the critical value 
Uf c the dynamics of the system is qualitatively different, 
reflecting the change in the behavior of the effective po- 
tential inversion points, see equation (|35[) . Let us start 
discussing the dynamics of double occupancy. Since the 
effective potential T (D) has two simple roots, the motion 
of double occupation D(t) is still periodic. However, the 
period T and the amplitude A of these strong coupling 
oscillations decrease upon increasing the strength of the 
quench, in contrast to the weak quench case. Indeed the 
latter simply reads A = Di — D + ~ 1/uf while the period 
reads 



T 



AK(l/k) 



with argument 1/fc given by 



2Di ufc 



(39) 



(40) 



k V u f i u f - u i) ' 

Deep in the strong coupling regime, Uf ^> u^, we get 

r^, (4i) 



smoothly matching the atomic limit result. Hence the 
resulting dynamics shows very fast oscillations with a 
reduced amplitude. In the strong quench limit the double 
occupation dynamics is completely frozen, doublons have 
no available elastic channel to decay^. 

As the critical quench line Uf c is approached from 
above the period of oscillations shows the same logarith- 
mic singularity found on the weak-coupling side. From 
equation (|39[) we immediately see that 



T 



log 



1 



Uf - Ufc 



(42) 



namely the same singularity, with the same prefactor, 
appears on the two side of the dynamical transition. 

As already anticipated, it is interesting to discuss the 
dynamics of the phase <p(t) when the quench is above the 
critical line. As shown in figure 01 as soon as the crit- 
ical line is crossed, the phase starts precessing around 
the whole circle (0, 2ir). This transition from a localized 
phase with small oscillations around <j) — to a delo- 
calized phase where the dynamics is unbounded is, from 
a mathematical point of view, completely analogous to 
what happen in a simple pendulum. Right at the critical 
quench line the dynamics is on the separatrix and the 
phase takes infinite time to reach its metastable config- 
uration. As we are going to see in the next paragraph 
this metastable configuration corresponds to a featureless 
Mott Insulator. Before concluding, let's briefly discuss 
the dynamics of quasiparticle weight Z(t) in the strong 
quench regime. As we see in figure similarly to the 
double occupation, also Z(t) shows fast oscillations with 
a period T given by (|4*T1) at strong coupling. Interest- 
ingly, the amplitude Az of those oscillations goes all the 
way to zero and keeps finite even for very large u / . This 
can be easily understood by looking at the dependence 
of the quasiparticle weight from the phase 4>{t). At half 
filling this simply reads (|27l) 



Z(t) = 16 D(t) (1/2 -D(t)) cos 2 <j){t), (43) 

from which we can conclude that, although the double 
occupation is neither zero nor one half, the quasiparti- 
cle weight can vanish due to its phase dependence. As 
a result of this vanishing minimum we conclude that for 
Uf ^> 1, even though the dynamics of double occupancy 
gets frozen in the initial state, the amplitude of oscilla- 
tions for Z remains constant and equal to Az = 1 — uf. 



Critical Line 

Quite interestingly, the weak and the strong coupling 
regimes that we have so far discussed are separated by 
a critical quench line Uf c at which mean-field dynamics 
exhibits exponential relaxation. This can be seen explic- 
itly since in this limit the effective potential is simply 
given by T (D) — D^j2uf c (Di — £>), thus the dynamics 
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FIG. 5: Dynamics after a quench at the critical interaction 
Uf c , for different initial conditions m. Both double occupa- 
tion D(t) and quasiparticle weight Z(t) decay exponentially 
to zero with a relaxation time r* ~ \j\fZi which increases 
with Ui approaching the initial Mott Insulator. In the bot- 
tom panel we compare the Gutzwiller results with those of 
DMFT (points, from Ref. 22 ) for a quench starting from the 
non interacting limit. 



is interesting to add some considerations. From our re- 
sults we see that for quenches at the critical line u f c the 
system reaches a steady state featuring a complete sup- 
pression of charge fluctuations, namely D — and Z = 0. 
This suggests that the above critical line Uf c is obtained 
by tuning the initial energy Eq of the quenched correlated 
metal to the energy of a collection of decoupled half filled 
sites, the ideal tij = Mott insulator. Indeed from this 
condition we immediately get 

1 + Ui 

Eq (Uf c , Ui) = E Mo tt — > u fc = — - — . (45) 

Surprisingly enough we find that the above condition 
gives a remarkable good agreement for the dynamical 
critical point found in DMFT. Indeed if we use that lat- 
ter criterium, we find an estimate for the critical Ut c in 
units of the hopping integral t and strating from Ui = 0: 

[/ /c = 4|£ fciTi |~3.3, (46) 

where E^in is the energy of a half-filled Fermi sea with 
a semielliptic density of states. Eq. (|4"6"|) is surprisingly 
close to the result of Refs. |22||23|. 



D. Long-time Averages 



can be easily integrated to obtain the double occupation 
D(t) at the critical quench line, 



D(t) = Di (1 - tanh 2 (t/r*)) . 



(44) 



We notice that, independently on the initial value of the 
correlation Ui, for Uf — Uf c the double occupancy re- 
laxes toward zero with a characteristic time scale r* = 
A/y/Zi that increases upon approaching the Mott insu- 
lator Ui — > 1. Analogously, also the quasiparticle weight 
Z(t) approaches zero for long-time, with the same expo- 
nential behavior, 

Z{t) = Zi (1 - tanh 2 (t/ T *)) . 

Since this is the only case in which our mean field dynam- 
ics features a long-time steady state it is worth to com- 
pare the above behavior to the DMFT result o 22 ' 23 . In 
figure [5] we plot the behavior of the quasiparticle residue 
Z{t) in the two approaches for the case Ui = 0. As we see 
they both vanishes at long times with a quite good agree- 
ment on the time scale. A similar comparison cannot be 
done for the double occupation D(t) which vanishes at 
long times in our mean field theory while saturates to a 
finite small value in DMFT. This is not surprising but 
again reflects the fact that our mean field dynamics can- 
not capture the role of incoherent excitations. The long 
time vanishing of th e q uasiparticle weight has been in- 
terpreted in Refs. 2'2 '2'A as a signature of thermalization. 
Although we cannot comment on this issue, since our 
mean field theory cannot account for thermalization, it 



As we have seen so far, the mean field Gutzwiller dy- 
namics is periodic in the main part of the phase diagram 
excluding the quench to the critical value ut c where an 
exponential behavior emerges. In spite of that, it is how- 
ever worth to investigate a properly defined long-time be- 
havior of the dynamics which, as we are going to show, 
features many interesting properties. To this extent we 
firstly introduce, for any given function 0(t) an integrated 
(average) dynamics defined through 

(0) t = - f dt'0{t') . (47) 
1 Jo 



Then it is natural to define the long-time average as 



O = lim t 



(0)t 



(48) 



Notice that, since the relevant observables are periodic 
functions of time with period To admitting a Fourier de- 
composition the above definition (|48[) can be equivalently 
written as 



0=^r [ dtO(t). 
lo J To 



(49) 



We now study the behavior of steady state averages as 
a function of the initial and final values of the interaction. 
We consider the half filled case and for simplicity we as- 
sume Uf > u t . Using equation (|49|) the average double 
occupation D can be written as which reads 



*4 



Di DdD 



(50) 



9 



where Di and D + has been defined in the previous sec- 
tion. In addition, due to energy conservation, the knowl- 
edge of the average double occupancy D completely fixes 
the average quasiparticle weight which reads 



Z = Zi + 8u f (D - D 



(51) 



We now evaluate the long-time average D and Z as 
given in Eq. 



(|50H51[) in the two different dynamical 
regimes we have previously identified. 



Weak Quenches: Uf < Uf c 

In the weak coupling regime and for Uf > itj the aver- 
age double occupation at long times reads 



D 



Di 1 



Ufc 
Uf 



D 



Uf c E(k) 



= Di 



uf K(k) 
Ufc (E{k)-K(k) 
u f 



V K{k) 



(52) 



where K(k) and E(k) are, respectively, the complete el- 
liptic integrals of the first and second kind with argument 



A: 



2 _ Ufjuf- 

2D: 



Similarly using the Eq. (fBTj) we get for 



the average quasiparticle weight the result 



Z = Zi 



m 

K(k) 



(53) 



It is interesting to consider the asymptotic regime of a 
small quantum quench Su — Uf — Ui — > 0. Then we can 
expand the elliptic integrals for small k to get 



D 



1 4 4 



(54) 



We see therefore that for small quenches the double oc- 
cupation follows the zero temperature equilibrium curve, 
independently on the initial value of the interaction Wj. 
This is clearly shown in figure [6] Since to lowest order in 
Su no heating effects arise, this result implies that after a 
small quench of the interaction the average double occu- 
pation D is thermalized. In addition, this result has an 
interesting consequence for what concerns the behavior 
of the quasiparticle weight Z. A simple calculation to 
lowest order in Su gives 



Zi -2uf (u f - , 



(55) 



from which we conclude that, as opposite to the double 
occupation l), the long-time average quasiparticle weight 
differs from the zero temperature equilibrium result even 
at lowest order in the quench Su. In particular if we 
evaluate Z for the special case of a quench from a non 
interacting Fermi Sea (ui = 0) for which Zi = 1 we get 
the result, 



I- Z(u f ) = 2(1- Z eq (u f )) 



(56) 




FIG. 6: Average double occupation D (top) and quasipar- 
ticle weight Z (bottom) as a function of uf at fixed Ui = 
0.0, 0.25, 0.5 compared to the zero temperature equilibrium 
result (dashed lines). 



firstly obtained in Rcf. 20 within the flow equation ap- 
proach. This peculiar mismatch between the zero tem- 
perature equilibrium quasiparticle residue and its non 
equilibrium counterpart is a general result of quench- 
ing a Fermi Se a 21 ' 33 . It signals the onset of a prether- 
mal regime where quasiparticle are well defined objects, 
momentum-averaged quantities such as kinetic and po- 
tential energy are thermalized while relaxation of distri- 
bution functions is delayed to later time scales. We note 
that our simple mean field theory correctly captures the 
onset of this long-lived state but fails in describing its 
subsequent relaxation towards equilibrium. 

Interestingly, when approaching the critical quench 
line from below the average double occupation D van- 
ishes logarithmically. Indeed for k — > 1 we have 

K{k) ~ log U/Vl - /c 2 ) + O (l - k 2 ) , (57) 



and 



therefore 



E (fc) ~ 1 + O (1 - k 2 



D~Di 



Uf-Uf c \ 2Di 



Uf 



loi 



\u fc — u f J 



(58) 



(59) 



The leading term is therefore logarithmic as mentioned, 
with linear corrections in Su = Uf c — Uf 



D 



2Di 



log ( -^—) 



Su log Su 



2u 



fc 



(60) 
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A similar behavior is found for the quasiparticle weight 
Z which reads 



2Z r 



log f -^—) 



(61) 



Strong Quenches: it/ > u/ c 

In the strong coupling regime the average double oc- 
cupation reads 



D = 



jfc-Uf Uf -Uj E(k) 
2 2 K{k) 



(62) 



with the argument given by k 2 — 2 f' Ufc . 

° ° J Uf{Uf—Ui) 

Deep in the strong coupling regime, Uf^>Ui,k goes to 
zero and we can use the asymptotic for E(k) and K{k) 



E(k) 
K(k) 



k 2 

1 ~ T 



to obtain 



D ~ Di 1 - 



2// 



(63) 



(64) 



We see therefore that, for infinitely large quenches, ut — » 
oo, the dynamics is trapped into the initial state. In- 
terestingly enough, for quenches starting from m = 
the scaling f|64[) exactly matches the strong coupling per- 
turbative result obtained in Ref. [22] for the prethermal 
plateau. Indeed using the fact that for m = we have 
DiUf c = 1/8 = |e|, where s is the kinetic energy of the 
half-filled Fermi Sea, we find 



D~Di 



21V 



in accordance with strong coupling perturbation the- 
ory. The agreement at strong coupling is remarkable 
if thought from the point of view of thermal equilib- 
rium, where one knows the Gutzwiller wavefunction can- 
not capture the Hubbard bands, and suggests that our 
Gutzwiller ansatz can interpolate between the weak and 
the strong coupling dynamical regime. 

As opposite, when approaching the critical quench line 
from above we obtain a vanishing long-time average, with 
the same logarithmic behavior we have found on the weak 
coupling side. Indeed for Uf —> Uf c from above we have 
that k —¥ 1~ and therefore we can again make use of the 
asymptotic for the complete elliptic integrals. We thus 
obtain 



D 



2Di 



\Uf-U fc ) 



lOf 



1 



8u log Su 
4A 



(65) 



Note that the approach to zero is the same in both sides 
of the phase diagram, while the corrections are slightly 
different. 



For what concerns the quasiparticle weight Z to get the 
leading behavior o(l/uf) we need the double occupancy 
to next-to-leading order. Expanding the ratio between 
elliptic functions we get 



E(k) 
K(k) 



h 2 k 4 

l-y-y + 0(fc 6 ), (66) 



and using the expression for k ~ Zi/Au 2 we obtain the 
following asymptotic behavior for Z 



2uj k 2 



l + fc 2 /4)-f 1 



Z, 



16u 2 



(67) 



which shows that also Z increases from the critical line to 
large u / and deep in the strong coupling regime it satu- 
rates to a finite plateau which, however does not coincide 
with its initial value Zi but rather it is smaller by a factor 
of two due to energy conservation. 



E. Quench Dynamics away from half-filling 

An important outcome of previous sections has been 
the identification of a critical interaction quench Uf c 
where an exponentially fast relaxation emerges. This 
value of quenches separates two different dynamical 
regimes where the system gets trapped into metastable 
prethermal states. In order to understand the origin of 
such a sharp transition and its possible relation to equi- 
librium critical point of the Hubbard model it is natural 
to extend the mean field analysis away from half-filling, 
where no transition between a Metal and a Mott Insula- 
tor exists in equilibrium. This can be done straightfor- 
wardly, for example, by a direct integration of the mean 
field equations of motion (|28H29|) . It is however more 
instructive to proceed again by considering the effective 
dynamics for the double occupation, obtained using the 
conservation of energy, that we wrote as 



D=y/T(D). 



(68) 



We now argue that any finite doping 8 is enough to 
wash out the dynamical critical point and turn it into 
a crossover. To see this, it is worth to consider again the 
effective potential T (D) which enters the above dynam- 
ics. Indeed the qualitative analysis we have performed 
in section Till CI can be done even for finite doping 5. As 
we will show explicitly in the Appendix [A] the effective 
potential keep the same structure for (5^0, with three 
inversion points respectively given by Di - the zero tem- 
perature finite doping Gutzwiller solution - and D±. As 
a consequence, all the differences between the doped and 
the half-filled case are hidden in the behavior of the two 
non-trivial roots Z3+,Z)_ as a function of Uf. Their ex- 
plicit expression is quite lengthy and it is reported for 
completeness in Appendix [XJ As we can see from figure 
[51 those two roots, which at half- filling are degenerate at 
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FIG. 7: Top Panel: amplitude Ad (left) and period Td at 
m = and S = 0.0,0.05,0.10,0.20,0.30. Bottom Panel: av- 
erages double occupancy D (left) and quasiparticle weight Z 
(right) at m = and 6 = 0.0, 0.05, 0.10, 0.20, 0.30. 



Ufa, are always distinct at finite doping. In particular at 
the half- filling critical quench line Uf c , we find at finite 
doping 

D + (uf c ) - D_(u/ C ) ~ S . 

As a consequence the dynamics of double occupancy (and 
hence of quasiparticle weight) always features a finite pe- 
riod given by 



~ y/u f (Di-D-) ' 

with the argument k of the elliptic function defined in 
term of the inversion points as 

k = V(A - D+) I (A - DJ) . (70) 

Notice that, since the two inversion points never collapse 
D + > D-, the argument k is always strictly lesser than 
one, k < 1, and no singularity in T arises. 

In figure [7] (top panels) we plot the period T and the 
amplitude A of the double occupancy oscillations in the 
doped case, as a function of ut at fixed it,. We notice 
that both quantities are smooth across it/ c , and in partic- 
ular the logarithmic singularity in the period turns into a 
sharp peak which broadens out as the doping increases. 

We finally remark that a small doping not only affects 
the dynamics, but also drastically changes the long-time 
averages properties with respect to the results we have 
depicted in section UlI CI This can be worked out explic- 
itly by using the same equations we have obtained for the 
half-filling case, cfr. section MI D[ provided the correct 
expression for the roots D +1 D- is used. As we can see 



from figure [7] both double occupation and quasiparticle 
weight stay always finite as «/ increases and only show 
a dip around the critical quench line which is gradually 
smoothed out as the doping increases. 

In conclusion, we have shown that the dynamical tran- 
sition described in section IIII CI is a peculiar feature of 
the half-filled case, namely that any finite doping S ^ 
is enough to wash out this dynamical transition, cutting 
off the logarithmic divergence in the oscillation period T. 



F. Discussion 

We conclude this section by discussing the results of 
our time dependent mean field theory for the fermionic 
Hubbard model in light of those recently obtained in the 
literature using different approaches, such as the Flow 
Equation method 20 i 21 and the Non Equilibrium Dynam- 
ical Mean Field Theor y 22 i 23 , both of which considered a 
quantum quench starting from an half-filled non interact- 
ing Fermi Sea. As we already mentioned, our mean field 
results feature an oversimplified periodical dynamics that 
lacks relaxation to a steady state at long times. This can 
be traced back to the suppression of quantum fluctua- 
tions which is at the ground of our treatment. In this 
respect we notice that both approaches work much bet- 
ter, displayng some damping at long times. Beside this 
obvious drawback we can say that, quite remarkably, a 
mean field theory catches many interesting features of 
the problem. 

First of all our variational ansatz is able to capture 
both regimes of pre-thermalization found at weak 20 and 
strong coupling 22 . Those long lived metastable regimes, 
which are, respectively, due to Fermi statistics and to 
long-lived double occupations, are quantitatively repro- 
duced by our approach as it appears clearly from the 
analysis of long time averages (see Eqs. ([55]) and (|M)) ). 
However, as generally expected in mean field theories, 
those metastable states are wrongly predicted to have 
infinite lifetime. A second interesting point that clearly 
emerges from our analysis is the existence of a dynamical 
critical line that separates those two distinct regimes, and 
where an exponentially fast relaxation emerges, as firstly 
shown in Ref. [22|. On one hand, the existence of a dy- 
namical critical point could be anticipated since at equi- 
librium the model undergoes a quantum phase transition, 
the Mott transition. Indeed, as we have shown, any finite 
doping turns the dynamical transition into a crossover. 
On the other hand, it was noted in Ref. 22 that the energy 
pumped in the quench at U f c with {/, = would corre- 
spond, should thermalization be assumed, to an effective 
temperature higher than the Mott ending point, where 
no critical dynamics could have been foreseen. Such an 
observation points to a dynamical transition that could 
be associated with loss of ergodicity and which is not in- 
compatible with our finding that the critical quench oc- 
curs when the correlated metal is initially prepared with 
the energy of the ideal Mott insulator, a collection of in- 
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dependent sites. Interestingly enough, such a condition 
gives a excellent match with the DMFT estimate of the 
dynamical critical point (see equation |46|) . Finally, we 
note that the issue of a non equilibrium dynamical tran- 
sition in the quench dynamics of interacting quantum 
systems seems to be of more general interest. Indeed re- 
cent investigations on the fully connected Bose Hubbard 
model 2 ^ and the scalar <p A mean field theory — reveals 
that a very similar phenomenon is present in these mod- 
els as well. Whether this is an artifact of the mean field 
approximation or rather a generic feature of the quench 
dynamics of interacting quantum systems in more than 
one dimension, as recent works would suggest^ is an in- 
teresting subject that requires further investigations. In 
this respect an interesting question is the role played by 
small quantum fluctuations on such a dynamical transi- 
tion. We will try to partially address this issue in the 
remaining part of this paper. 



IV. Z 2 SLAVE SPIN FORMULATION 



A minimal formulation in terms of a single Ising spin 
and a fermionic degrees of freedom has been recently 
introduced , 25 ! 26 based on a mapping between the local 
physical Hilbert space of the Hubbard Model and the 
Hilbcrt space of the auxiliary model subjected to a con- 
straint. Here we derive the same mapping by showing 
that the identification holds for the partition functions 
as well, when evaluated order by order in perturbation 
theory in U . The advantage of this alternative formula- 
tion is that the role of the lattice coordination emerges 
more clearly. 

We write the Hubbard interaction as 



u r 



2(n-l) 2 



+ j(2n-l). 



The last term can be absorbed into the chemical poten- 
tial, so that we shall consider as interaction only the first 
term. We define 



We have shown that, within the Gutzwiller approx- 
imation, the variational principle when applied to the 
Shrcedinger equation amounts to determine the saddle 
point of an action S [4>i a , Oi a , $] that depends on pairs 
of conjugate fields, <fii a and Oi a , and on a Slater deter- 
minant or BCS wavefunction. The saddle point reduces 
to a set of first order coupled differential equations for 
the conjugate fields and for the average values of single 
particle operators on One could be tempted to 

interpret this result as the mean field decoupling of the 
Heisenberg equations of motion for the average values of 
a set of quantum operators corresponding to some effec- 
tive quantum Hamiltonian. Identifying such a quantum 
Hamiltonian could then allow adding quantum fluctua- 
tions on top of the mean field results. This is right the 
same conceptual scheme invoked to associate the time- 
dependent Hartree-Fock equations to an effective Hamil- 
tonian of non-interacting bosons that represent particle- 
hole excitations. In our case we would expect the quan- 
tum Hamiltonian to describe free electrons coupled to 
a set of conjugate Bose fields, 4n a and Oi a , which in 
fact resembles the conventional slave-boson approaches 
to correlated systems. 

We are going to show that this program can be easily 
accomplished in the simple Hubbard model, although in 
a different and more rigorous manner than simply quan- 
tizing the classical equations of motion. To this extent 
we formulate the original Hubbard model in terms of an 
auxiliary Quantum Ising Model in a transverse field cou- 
pled to free fermionic quasiparticles, in the framework of 
the recently introduced Z2 slave spin theor y 25 ' 26 . 

A. Mapping onto a Quantum Ising Model in a 
Transverse Field 

The idea of writing the Hubbard model in terms of aux- 
iliary spins coupled to free quasiparticles is not new i 34 ' 35 



2(n- l) 2 



n, 



(71) 



where the operator Q is real and unitary and has eigen- 
values — 1 for n = 1 and +1 for n = 0,2. It follows 
that 



(72) 



namely it changes sign to the fermion operator. 

Let us concentrate on a given site, with local energy e, 
whose Fermi operator we shall denote as c' a and density 
operator n. The rest of the lattice sites, Fermi oper- 
ators dR CT , are described by the generically interacting 
Hamiltonian Hbath and are coupled to the site under in- 
vestigation by 



H, 



We shall denote as 



■E 

Rcr 



*R c U Ra + H.c. 



(73) 



Hn — Hi 



bath 



Ht; 



the unperturbed Hamiltonian and 



4 4 ' 



the perturbation. Suppose we calculate the partition 
function within perturbation theory. A generic n-th or- 
der correction to the partition function is 
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z "" - (-?)" [ * r *> 



dr n Tr 



e -(/9-ri)i?o ^ e -(n-r 2 )J?o ^ ft e -(r„_i-r„) H Q e ~(r„-0) H 
I 



Because of (J72J| fJ i?o ^ = #6at/i + e n — H tun n = Hi . We even case one easily realizes that 
shall distinguish the two cases of n even or odd. In the 

I 



Z 



(2n) _ 



f7 



2n ,.0 



n JO 



dr 2n 



Tr 



e -(/3-n)J? e -(ri-T 2 )Hi e -(T 2 -r 3 )i? e -(r 2 „-i-T 2 „) Hi 6 -(t 2 „-0) H 



(74) 



r 



which resembles an iterated X-ray edge problem, like in We note that, since fi = 1, Eq. (|74[) is invariant under 
the Anderson- Yuval representation of the Kondo model. Hq f> Hi. For the odd case, one finds instead 

I 



Tr 



e -(/8-T"i) H g-(Ti-T 2 ) Hi e -(r 2 -T 3 ) H g-(T 2 „_i-T 2n ) Hi e -(T 2n -r 2n+ i) H e -(T 2n+1 -0) Hi ^ 



(75) 



Once again the above expression is also equal to that one 
where H is interchanged with Hi . 

Can one reproduce the same perturbative expansion 
with some other model? Let us consider an Ising-like 
Hamiltonian Hi S i ng = H * + H tra nsv where the unper- 
turbed term is 



H* = H bath + en + a x H t , 



the perturbation is 



(76) 



(77) 



and a a , a = x,y,z, are Pauli matrices. If we take the 
trace over eigenstates of a x - note that for a x = 1 H* = 
H , while for a x = — 1 i?» = Hi - the perturbation ([77j) 
may act only an even number of times and one easily find 
that the final result is just twice (|74p . In other words, 
2;(2n) - 1S ^gjf £ 2n-th order term in the perturbative 
expansion of the Ising model Hj sing . How do we get the 
odd order terms in the expansion? Let us consider the 



perturbative expansion of 



-Tr e-^ Hlsi "o a z fl 



It is clear that now only odd terms in the expansion over 
eigenstates of a x will contribute and one easily realizes 
that the final result is twice (I75|) . 

Therefore, the partition function of the original model 
is also equal to 



Z = Tr 



-pu ls 



i-a z n 



We note that 



Q = 



(78) 



(79) 



is actually a projector of the enlarged Hilbert space onto 
the subspace where if n = 1 then a z = + 1 while, if 
n = 0,2, then a z = — 1. As a matter of fact, Q is just 
the constraint introduced in Ref. 26i as a basis of the Z2 
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slave-spin representation of the Hubbard model. In fact, 
what we have done here is simply re-deriving the mapping 
of Ref. in a different way. There are however some 
interesting aspects of the mapping that emerge clearly at 
the level of the partition functions and were not discussed 
in Ref. SI 

We note that what we have shown so far is that, given 
an Anderson impurity model with Hamiltonian 



IM 



Hbath + Hfann +£ll+y(n-l) 2 

= H bath + H tunn + e n + — (1 + fi) , (80) 
its partition function can be also written as 



Z A 



1 



IM 



Tr 



- Zi s i n g ( 1 - (<J Z CI) 



(81) 



where 



U 



Hi sing = H ba th + £ Tl + CT X H tu nn + ~£ (1 - CT Z ) , (82) 



and 



Zisina = Tr ( e -^« 



As mentioned above, Zi S i ng is even in U, while (a z fi) is 
odd. As a simple byproduct, we note that, if particle- hole 
symmetry holds, the partition function must be even in 
U, so that 



Z a 



IM 



I sing i 



(83) 



hence the constraint is uneffective and the mapping holds 
trivially. It was noticed in Ref. 26 that H] S i ng in l|82p 
possesses a local Z2 gauge symmetry, cj — > — cj and 
a x — > —<t x , which can not be broken. Indeed, the factor 
1/2 in (|8"3"|) avoids the consequent double counting. 

One can straightforwardly extend the above procedure 
to a collection of interacting sites, hence to the Hubbard 
model, with the final result that 



Tr 



Se 



-PBt 



(84) 



where now 

Hi sing = -t C7 R (T R' c L c rv + jE (!-or)- 

R 

(85) 



<R,R'> a 

and the constraint is 

R ^ ' 

We note that, if t — > —it, the mapping still holds and 
shows that the time-evolution of the Hubbard model can 
be mapped onto the time evolution of Hj sing . In partic- 
ular, since [Q, Hj sing ] = 0, the two evolutions are exactly 
the same on a state that satisfies the constraint. 



B. Recovering the Gutzwiller approximation at 
equilibrium 

Let us now consider a lattice whose coordination tends 
to infinity in a such a way that the hopping energy per 
site remains well defined. In this limit, it is well known 36 
that the Hubbard model maps onto an Anderson im- 
purity model self-consistently coupled to a conduction 
bath. We showed earlier that when particle-hole symme- 
try holds, the constraint is uneffective for the mapping 
of the Anderson impurity model to the Ising model. It 
follows that the same holds also for the Hubbard model, 
in which case 



Z 



Hubbard 



N 



Zls 



(87) 



where N is the number of sites. 

Therefore, in infinite coordination lattices and at 
particle-hole symmetry, we could calculate the partition 
function of the model 



H 



t 



Ising 



u 



<R,R'> a 



R 



and obtain that of the Hubbard model through (f87| . The 
factor z in (|88|) is the lattice coordination and must be 
sent to infinity at the end of the calculation^ 6 - It turns 
out that the Gutzwiller approximation is nothing but the 
mean field decoupling of Hj S i ng , assuming a wavefunction 
product of an Ising part times a fermionic one. The de- 
generacy of the solution that derives from the local Zi 
gauge symmetry, crf^ — > — cf^ and cj^ CT — > — c^ CT is can- 
celed out by the (1/2) N factor in (|57)l. 

To recover the Gutzwiller result for the Mott transi- 
tion, let us consider a trial translationally-invariant wave- 
function \^>) — |$ CT ) |$ c ), where is a Ising-spin state 
and |$ r ) an electron one. If we define 



3 Rcr C R'(T 



H.C.\$ C ) = — 



where —e is the average hopping energy per site of 1$ 
then the average value per site of the Hamiltonian 
is 



E=(®„\ -e 



- E 

<R,R'> 



R 



i.e. the energy of an Ising model in a transverse field. 
We assume = U |$o) where the unitary operator 



w = exp (4 E CT ^) 



(89) 



so that E becomes the average value on |$o) of the Hamil- 
tonian 

H* = ^ 1 ~ cos /3 a^- sin ^a^ 



R 
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e- E (cosjScT-R-sinjSo^,,^ 

> 

(cos/3er R , - sin/3o- R ,). (90) 



<R,R'> 



the -Z2-slave-spin theory for correlated fermions, to which 
we refer for further details. In what follows, we shall in- 
stead discuss a way to add quantum fluctuations in an 
out-of-cquilibrium situation. 



We assume that |$g) is so close to the fully ferromagnetic 
state with all spins oriented along x that we can set 



"ft 



1- (4 



i - n 



R- 



(91) 
(92) 
(93) 



where ir and p R are conjugate variables. If we substi- 
tute the above expressions in (j9"0)) and fix j3 in such a way 
that all terms linear in ir vanish, we find 



• R U 

8e 



(94) 



for C7 < 8e = C/ c , while sin/? = 1 otherwise. t/ c is the 
mean-field value of critical transverse field that separates 
the ordered phase from the disordered one in the Ising 
model. It also identifies the Mott transition in the origi- 
nal Hubbard model, and, in fact, the value of U c coincides 
with that of the Gutzwiller approximation. Because of 
the above choice of /3, once we expand the Hamiltonian 
(1901) up to second order in cc R and p R we find, apart from 
constant terms and in units of U c , 



a x ■> 
2 4^ 



(*R 



PR 



b 2 

2 z 



<R,R'> 



(95) 



where a = 1/2 and 6 = u 2 /2 for u = (U/U c ) < 1, the 
metallic phase, while a — u/2 and b = 1/2 for u > 1, 
the Mott insulator. The spectrum of the excitations on 
both side of the transition is that of acoustic modes with 
dispersion in momentum space 



w q = \ a (a - 67, 



(96) 



where, assuming a hypercubic lattice in d — z/2 dimen- 
sions, 



1 

7q = ^E cosqa e t -1 ' 1 ] 



(97) 



with q a the components of the wavevector q. At the 
transition a — b and the spectrum becomes gapless at 
q = 0. In principle, at the same level of approximation 
one should also take into account the coupling between 
the spin-waves of the Ising model and the conduction 
electrons via the hopping term in (|88|) . We just mention 
that, deep in the insulating side, where u; q ~ u/2 ^> 1, 
one can integrate out the acoustic modes and obtain 
the antiferromagnetic Heisenberg model known to be the 
large U limit of the half-filled Hubbard model. A thor- 
ough analysis of the role of quantum fluctuations at equi- 
librium has been presented in Ref. [26j in connection with 



C. Recovering the Gutzwiller approximation 
out-of-equilibrium 

Because the two models can be mapped onto each 
other, a quantum quench in the Hubbard model is equiv- 
alent to suddenly change the transverse field in the Ising- 
like model ([55)1 at particle-hole symmetry and in the 
limit of infinite coordination lattices. We shall keep as- 
suming a factorized time-dependent trial wavefunction 
|$ CT (f)) |$ c (*)), each component |$o-(f)) and |$ c (*)) be- 
ing translationally invariant. The electron wavefunction 
will evolve under the action of a time-dependent hopping, 
which is however still translationally invariant. Hence, if 
|$ c (t = 0)) is eigenstate of the hopping at t < 0, in 
particular its ground state state, it will stay unchanged 
under the time evolution. Therefore we shall only focus 
on the evolution of the Ising component. Its Hamiltonian 
at positive times and in units of U c is 



R 



1 2 x ^ 

<RR'> 



CT R CT R' 



(98) 



and we assume that at time t — \$> a (t = 0)) is the 
approximate ground state defined in the previous sec- 
tion IIV Bl for a different transverse field Uj . The time- 
evolution is thus described by the Schrcedinger equation 



We assume 



where now 



id t \^(t)) = H\^(t)). 



|<M*)>=w(*)l*o(*)>, 



(99) 



(100) 



U(t) = exp i 



a(t) 



(101) 

It follows that |$o) must satisfy the equation of motion 

i3 t I $0 (*)}=#*(*) I $o(*)), (102) 
where, apart from constants, 

H*(t) = -iU{tfU{t)+U{tf HU{t) (103) 



E 



R L 



- cos^cr R -- sin P CT R + -cr R 



cos a cos P cr R + cos a sin p cr R 



sin a a 



R 
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^ 2 / 

g - J! (cos£cr R -sin/3cr R 



^ cos [3 cr R , — sin /3 <t r , 



<R,R'> 



be dangerous, as we shall see. We shall analyze the time 
dependent problem (|106[) separately in the two different 
cases of quenching from the correlated metal or from the 
Mott insulator, starting from the latter that is simpler. 



In the same spirit of the spin-wave approximation above, 
we shall assume that |$o(^)) is a t any time close to a 
fully polarized state along x, so that we can safely use the 
approximate expressions (|9Tj) - (|93l) for the spin operators. 
Just like before, we fix a(t) and j3(t) in such a way that all 
linear terms in xr and pr vanish and find the following 
set of equations 



p = — ^- ana 



f U f 

- cos (3 — r- cos a cot (3. 



(104) 
(105) 



These equations have to be solved starting from the ini- 
tial condition appropriate to the approximate ground 
state with transverse field Ui, i.e. a(0) — and sin/3(0) = 
Ui if Ui < 1 otherwise sin /3(0) = 1, see Eq. (J94j . In ad- 
dition, as noticed before, the equations admit a constant 
of motion, which can be regarded as the classical energy, 



E 



u f o 1 2 a 

- — f- cos asm/)-- cos p. 
A 8 



One can readily recognize that the dynamical system 
(|104p - (|105p is equivalent to that one we previously ob- 
tained within the time-dependent Gutzwiller approxima- 
tion. However, as we are going to see in the next section, 
this alternative formulation however allows us to access 
quantum fluctuations, assuming they are small. 



D. Quantum Fluctuations beyond mean field 
dynamics 

The time dependent hamiltonian i?*(t) we have ob- 
tained in the previous section, Eq (|103[) . accounts in prin- 
ciple for quantum fluctuation effects. A simple way to 
proceed is to fix the parameters a(t) and j3(t) in such a 
way Eqs. (|104|) and ([105JI arc satisfied, and expand the 
hamiltonian up to second order in xr and pr. The result 
has no more linear terms and simply describes coupled 
harmonic oscillators with time-dependent parameters. 



#,(*) 



Uf cos a(t) 
A sin /3 (t) 



R 



sin 2 (3(t) 2 
A z 



E 

<R,R'> 



(106) 



We note that such a treatment, similar to what we have 
done in equilibrium, is equivalent to include gaussian fluc- 
tuations without renormalizing the transition point. In 
other words, we are studying the effect of quantum fluc- 
tuations around the semiclassical trajectory without al- 
lowing any feedback of these on the latter, which could 



1. Quenching from the Mott insulator 

In this case Ui > 1 and the initial values of the Euler 
angles are a(0) = and sin/3(0) = 1. It follows from 
Eqs. (|105[) and (|104[) that these angles will not evolve in 
time so that H* in (I106P does not depend on time and co- 
incides with for a — ut/2 and b = 1/2. This Hamil- 
tonian is well defined provided Uf > 1, which simply 
reflects that our assumption of weak quantum fluctua- 
tions loses its validity if the quench is too big. Therefore 
we shall assume Uf > 1, namely a quench withing the 
Mott insulator domain. 

Initially the system is described by the Hamiltonian 
with a — Ui/2. We assume that the initial state is 
the ground state of such a Hamiltonian. At times t > 0, 
this state is let evolve with the same Hamiltonian, but 
now with a = Uf /2. This problem can be readily solved, 
being equivalent to starting from the ground state of a 
harmonic oscillator and evolving it with a Hamiltonian 
having different mass and spring constant. We find that 
the time-dependent average value of the double occu- 
pancy is 



D(t) 



16V" 



1 



K K 



K 2 



>q 



K 2 

Kin 



K 



COS 2iOqt 



,(107) 



where 



see (|96J and ((971), while 

Ui 



K 



Ui - 7q 



'/ 



u f -7q 



are the parameters of the canonical transformation to 
find the normal modes of the initial and final Hamilto- 
nians, i.e. x — >■ y/K x and p — > p/yK. Seemingly, the 
hopping renormalization factor Z(t) turns out to be 

Z(t) - (ofof) 



2V^2 7q 
q 




Ki, 



K 



/q 



Ki, 



COS 2ujqt 



(108) 



We note that the sum of the oscillatory terms in (|107[) 
and (|108[) vanishes for t — >• oo, unless ut — > oo, so that 
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tU u. 

C 1 



FIG. 8: Behavior of the frequency function of time 

for Ui =0.1 and Uf = 0.2 (top panel) and u/ — 0.6 (bottom 
panel). We see that for suitable values of Uf the frequency 
can become negative for some time intervals. 



asymptotically D(t —> oo) and Z(t —> oo) approach val- 
ues that do not corresponds either to the initial ones nor 
to the equilibrium values for u — Uf. 

We remark that the above time evolution derives just 
by the quantum fluctuations. Should we neglect these 
latter, we would not find any dynamics for these quanti- 
ties. 



2. Quenching from the metal 

We now consider the case in which Ui < 1 so that 
initially a(0) = and sin/3(0) = Uj. With such initial 
values, the time evolution controlled by (I105P and (|104[) 
is non trivial, unlike the previous example of a Mott insu- 
lating initial state. As we mentioned before, the Hamilto- 
nian describes coupled harmonic oscillators with time 
dependent parameters. The time dependent frequency of 
these oscillations reads 



Uf cosa(i) 



sin 



7q 



(109) 



with 7 q defined in Eq. (|97p . Since the minimum fre- 
quency is obtained for q — we immediately realize 
that in order to have stable fluctuations the condition 
UfCosa(t) > cos 3 /3(t) has to hold. 

In figure [8] we plot the behavior of Lu^ =0 (t) as obtained 
from the semiclassical dynamics. We notice that for suit- 
able values of Uf it exist multiple time intervals at which 
w q=oM < an d fluctuations become unstable. In par- 
ticular, by looking at the mean field dynamics, it is easy 
to realize that there is a whole region of quenches, just 
around the dynamical transition, for which an instability 
in the fluctuation spectrum may occur. This region of 



FIG. 9: Behavior of the instability lines uju defined in the 
main text, as a function of < m < 1. We see that these 
lines bound a region of the phase diagram around the mean 
field critical line Uf c where fluctuations grows exponentially 
in time. This region shrinks upon approaching m — > 1 while 
becoming wider and wider in the opposite regime of quenches 
from a non interacting Fermi system. 



unstable modes is bounded by two lines u*^, u*^ 2 whose 
behavior is plotted in figure [9] The line u*f 2 can be ob- 
tained analytically by simple means and reads 



'/2 



2ui 



(110) 



As a result of this analyis we conclude that for quenches 
below and above these instability lines we can use the spin 
wave approximation to compute corrections to quantum 
dynamics, since all q-modes are stable. As opposite for 
quenches around the critical mean field line part of the 
spectrum becomes unstable. Conversely, we previously 
found that the same method is, at least, well defined 
when quenching from the Mott insulator down to the 
Mott transition. We believe that this difference is not ac- 
cidental and that the dynamics of quantum fluctuations 
quenching from the metallic side is poorly described by 
the Hamiltonian (|106[) . The metallic phase corresponds 
in our language to the ordered phase of the Ising model, 
where a finite order parameter is spontaneously gener- 
ated. The equations of motion (|104l) and (jl05[) describe 
the dynamics of the condensate alone. The approach 
in section TlV CI implicitly assumes quantum fluctuations 
that follow adiabatically the evolution of the condensate. 
However, the quantum fluctuations must in turn affect 
the evolution of the condensate, a feedback that is ab- 
sent in the above scheme and explains why the latter fails 
if the quench is big enough. Anyway, the fact that the 
Hamiltonian (|106[) become unstable before the dynami- 
cal critical point is encountered suggests that the effect 
of quantum fluctuations grows and it is not unlikely to 
modify substantially the dynamics. 
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V. CONCLUSIONS 

We have introduced a variational approach to strongly 
correlated electrons out of equilibrium. The idea is to 
give an ansatz on the time dependent many-body wave 
function and to obtain dynamical equations for the pa- 
rameters by imposing a saddle point on the real-time 
action. While this strategy is widely used for non in- 
teracting fermionic systems, in the spirit of time depen- 
dent Hartree-Fock, its extension to strongly correlated 
electrons represents a novelty with many possibilities for 
further developments. Applications of this method can 
range from dynamics in closed quantum systems to non 
equilibrium transport in correlated quantum dots, for 
which a related variational approach for the steady has 
been recently proposed^!. 

In this paper we have applied this variational scheme 
to the single band Hubbard model using a proper gen- 
eralization of the Gutzwiller wavefunction. It is worth 
mentioning, however, that the method is general and can 
be applied also to other correlated wavefunctions, as long 
as a suitable numerical or analytical approach is available 
to calculate the variational energy controlling the clas- 
sical dynamics of the variational parameters. As a first 
application we have studied the dynamics of the Hub- 
bard model after a quantum quench of the interaction. 
This is an interesting open problem for which results 
have been obtained only very recently using sophisticated 
non equilibrium many body techniques. Remarkably, al- 
though extremely simple, our approach seems to capture 
many non trivial effects of the problem and shows a good 
overall agreement with the picture provided by DMFT. 
From this perspective it can be seen as a simple and intu- 
itive mean field theory for quench dynamics in interacting 
Fermi systems. 
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Appendix A: Details on the Gutzwiller calculations 
at finite doping 



In this appendix we describe in some detail the analyis 
of the mean field dynamics at finite doping. We start 
from the equation (|33p where the phase <fi is expressed in 



terms of D using energy conservation 

E -UD- 2e{n - 2D) (y/W 



M{n-2D) (D + 5) 

(Al) 

This result can be inserted into the equation for D{t), 
which reads after simple differentiation 

= -8s (n - 2D) ^D (D + S) siri(f> coscf> . (A2) 

After some simple algebra we end up with a differen- 
tial equation for the time-dependent double occupation 
whose general structure is 



dD 
~~dt 



±vnDj, 



(A3) 



where T(D) can be thought as an effective potential con- 
trolling the dynamics of D(t). Its explicit expression 
reads T (D) = T + (D) T_ (D) where 

r±(D) = ± Eq-UD -2e(n-2D) (yiTTS±v r D 

After some lengthy but straightforward calculations it is 
possible to bring the function T (D) to a polinomial form, 
namely to 



r (£>) = 73 D 3 + 72 D 2 + 7i D + 7 o , 



(A5) 



where j a 's are coefficients depending on the initial Ui and 
final U f interactions as well as on the doping S. We first 
notice that for 5 — the expression for T simplifies to 
read 

T s=0 (D) = {u f D- E ) (E -u f D + 2D (1/2 - D)) , 

(A6) 

where the initial energy E$ reads as in Eq. (|32j) . It is 
easy to very that the effective potential has three roots 
Di, D±, the former corresponding to the equilibrium 
Gutzwiller solution at T = 0, Di = (1 — m) /4 while 
the latters given respectively by 



and 



D + = 



Uf < Uf c 



u f > u fc Di (l - 



u } < u fc Di 

. Uf c —Uf 

Uf > u fc 2 . 



In the doped case we cannot obtain expressions as simple. 
However we notice that T+ {Di) = 0, since by construc- 
tion 

E =UfDi+2e{n-2Di) UfD~+5 + v 7 ^) ■ (A7) 
As a consequence we can write the effective potential as 



r (D) = {D- Di) $ (£>) 



(A8) 
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with $(D) that can be formally written as 

$ (£>) = 73 D 2 + ( 72 + Am) D + (ji + D Ql2 + Dfe) 

(A9) 

once the definition of the effective potential as a polyno- 
mial in D, Eq (| A5|) . is considered. From this result we 
obtain for the other two inversion points D± the follow- 
ing result 



(72 + A7s)TV / A 



4u 



f 



(A10) 



with A = ( 72 + A73) 2 - 473 (71 + A72 + Df lz ). The 
explicit expression for the coefficients 7 a can be easily 
found after some simple but lengthy algebra. These read 



73 
72 

7i 



-2u 



f 



2E 



2u f E 



u f {l-26)-S 2 /4 
- £b (1 - 26) 



+ 



with Eq given by Eq. (IA7|) . It is interesting to note that 
all the dependence from the initial interaction is hid- 
den into the Gutzwiller equilibrium solution D. L . The 
qualitative analysis can proceed along the same lines as 
in the previous section, the only difference being that 
Di is not known analytically. By solving the equilib- 
rium Gutzwiller problem at finite doping (see appendix) 
we can easily obtain Dj, hence D T through Eq (|A10[) . 
When inserted back into the previous results for T and 
A we find further evidence that no singularity emerges 
for any finite S in those quantities, which nevertheless 
features some signature of the zero doping criticality. In 
particular both T and A are smooth functions displayng 
a sharp peak around Uf c . 
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